module smooth_indicators_mod

  implicit none

  interface
    pure real(8) function smooth_indicator_interface(a)
      real(8), intent(in) :: a(:)
    end function
  end interface

contains
pure real(8) function smooth_indicator_2(a) result(res)

  real(8), intent(in) :: a(:) ! 2

  res = a(2) * a(2)

end function smooth_indicator_2

  pure real(8) function smooth_indicator_3(a) result(res)

    real(8), intent(in) :: a(:) ! 3

    res = a(2) * a(2) + 13 * a(3) * a(3) / 3

  end function smooth_indicator_3

  pure real(8) function smooth_indicator_rhombus_2(a) result(res)

    real(8), intent(in) :: a(:) ! 3

    res = a(2) * a(2) + a(3) * a(3)

  end function smooth_indicator_rhombus_2

  pure real(8) function smooth_indicator_2x2(a) result(res)

    real(8), intent(in) :: a(:) ! 4

    res = a(2)**2 + a(3)**2 + 7 * a(4)**2 / 6

  end function smooth_indicator_2x2

  pure real(8) function smooth_indicator_3x3(a) result(res)

    real(8), intent(in) :: a(:) ! 9

    res = (  720 * a(2) * a(2) + 3120 * a(3) * a(3) + 720  * a(4) * a(4) + 840   * a(5) * a(5) &
          +  120 * a(4) * a(6) + 3389 * a(6) * a(6) + 3120 * a(7) * a(7) + 120   * a(2) * a(8) &
          + 3389 * a(8) * a(8) + 520  * a(3) * a(9) + 520  * a(7) * a(9) + 13598 * a(9) * a(9) ) / 720

  end function smooth_indicator_3x3

  pure real(8) function smooth_indicator_4x4(a) result(res)

    real(8), intent(in) :: a(:) ! 16

    res = a(2)**2 + (13*a(3)**2)/3 + (3129*a(4)**2)/80 + a(5)**2 + (7*a(6)**2)/6 + a(5)*a(7)/6 + (3389*a(7)**2)/720        &
          + (17./30)*a(6)*a(8) + (47459*a(8)**2)/1120 + (13*a(9)**2)/3 + a(4)*a(10)/24 + (3389*a(10)**2)/720               &
          + (13./18)*a(3)*a(11) + (13./18)*a(9)*a(11) + (6799*a(11)**2)/360 + (1043./160)*a(4)*a(12) + (73./32)*a(10)*a(12)&
          + (22846129*a(12)**2)/134400 + a(2)*(12*a(4) + 4*a(10)/24 + a(12)) + 0.5*a(5)*a(13) + a(7)*a(13)/24              &
          + (3129*a(13)**2)/80 + (17./30)*a(6)*a(14) + (11./80)*a(8)*a(14) + (47459*a(14)**2)/1120 + a(5)*a(15)/24         &
          + (73./32)*a(7)*a(15) + (1043./160)*a(13)*a(15) + (22846129*a(15)**2)/134400 + (11./80)*a(6)*a(16)               &
          + (114997*a(8)*a(16))/5600 + (114997*a(14)*a(16))/5600 + (19583517*a(16)**2)/12800
    
  end function smooth_indicator_4x4
  
  pure real(8) function smooth_indicator_5x5(a) result(res)
    real(8), intent(in) :: a(:)
    
    res = a(2)**2. + (13.*a(3)**2)/3. + (3129.*a(4)**2)/80. + (87617.*a(5)**2)/140. + a(6)**2. + (7.*a(7)**2)/6. + (1./6)*a(6)*a(8) + (3389.*a(8)**2)/720. + (17./30)*a(7)*a(9) + (47459.*a(9)**2)/1120. + (1./40)*a(6)*a(10) +    &
          (5101.*a(8)*a(10))/1120. + (54673043.*a(10)**2)/80640. + (13.*a(11)**2)/3. + (1./24)*a(4)*a(12) + (3389.*a(12)**2)/720. + (7./20)*a(5)*a(13) + (13./18)*a(11)*a(13) + (6799.*a(13)**2)/360. + (1043./160)*a(4)*a(14) +   &
          (73./32)*a(12)*a(14) + (22846129.*a(14)**2)/134400. + (87617./840)*a(5)*a(15) + (13./120)*a(11)*a(15) + (306967.*a(13)*a(15))/16800. + (469977913.*a(15)**2)/172800. + (1./2)*a(6)*a(16) + (1./24)*a(8)*a(16) +          &
          (1./160)*a(10)*a(16) + (3129.*a(16)**2)/80. + (17./30)*a(7)*a(17) + (11./80)*a(9)*a(17) + (47459.*a(17)**2)/1120. + (1./24)*a(6)*a(18) + (73./32)*a(8)*a(18) + (24721.*a(10)*a(18))/22400. +                             &
          (1043./160)*a(16)*a(18) + (22846129.*a(18)**2)/134400. + (11./80)*a(7)*a(19) + (114997.*a(9)*a(19))/5600. + (114997.*a(17)*a(19))/5600. + (19583517.*a(19)**2)/12800. + (1./160)*a(6)*a(20) +                            &
          (24721.*a(8)*a(20))/22400. + (37850569.*a(10)*a(20))/115200. + (3129.*a(16)*a(20))/3200. + (2105043.*a(18)*a(20))/12800. + (368483712607.*a(20)**2)/15052800. + (21./5)*a(11)*a(21) + (7./20)*a(13)*a(21) +              &
          (21./400)*a(15)*a(21) + (87617.*a(21)**2)/140. + (1./160)*a(4)*a(22) + (5101.*a(12)*a(22))/1120. + (24721.*a(14)*a(22))/22400. + (54673043.*a(22)**2)/80640. + (21./400)*a(5)*a(23) + (7./20)*a(11)*a(23) +              &
          (306967.*a(13)*a(23))/16800. + (7071./800)*a(15)*a(23) + (87617./840)*a(21)*a(23) + (469977913.*a(23)**2)/172800. + (3129.*a(4)*a(24))/3200. + (24721.*a(12)*a(24))/22400. + (2105043.*a(14)*a(24))/12800. +             &
          (37850569.*a(22)*a(24))/115200. + (368483712607.*a(24)**2)/15052800. + (1./480)*a(2)*(240.*a(4) + 80.*a(12) + 20.*a(14) + 12.*a(22) + 3.*a(24)) + (87617.*a(5)*a(25))/5600. + (21./400)*a(11)*a(25) +                    &
          (7071./800)*a(13)*a(25) + (2475532433.*a(15)*a(25))/940800. + (87617.*a(21)*a(25))/5600. + (2475532433.*a(23)*a(25))/940800. + (2210903809027.*a(25)**2)/5644800. +                                                      &
          (a(3)*(15120.*a(5) + 2600.*a(13) + 1260.*a(15) + 390.*a(23) + 189.*a(25)))/3600
                          
  end function smooth_indicator_5x5

end module smooth_indicators_mod
